rm(list = ls())

if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, rdrobust, pbapply)

# Source helper function to tidy rdrobust results

source("code/helper_function.R")

## Load data

forsa <- read_rds("data/forsa.rds") %>%
    filter(east == 1) %>%
    filter(!is.na(vote_afd)) %>%
    filter(!is.na(yob)) %>%
    filter(yob > 1945) %>%
    filter(year > 2016) %>%
    mutate(hh_income = as.numeric(hh_income))

glimpse(forsa)

## List of outcomes

outcomevars <- c(
    "hh_income",
    "unemployed"
)

labels <- c(
    "Household income", "Unemployed"
)

## If continuous, standardize

forsa <- forsa %>%
    mutate(hh_income = hh_income / sd(hh_income, na.rm = T))

## Make samples

ss_list <- list(
    forsa %>% filter(male == 1),
    forsa %>% filter(male == 0)
)
ss_labs <- c(
    "Men",
    "Women"
)

##

cutoffs <- c(1971:1972)

# Run models

res_main_m <- pblapply(1:length(ss_list), function(i) {
    lapply(cutoffs, function(cutoff) {
        lapply(outcomevars, function(o) {
            print(o)

            df_use <- ss_list[[i]]
            df_use[, "y"] <- df_use[, o]
            df_use <- df_use %>%
                filter(!is.na(y))

            ## W/o covs

            m1 <- rdrobust(
                y = df_use$y,
                x = df_use$yob,
                c = cutoff,
                deriv = 1,
                scalepar = 1,
                p = 1
            )

            tidy_rd(m1, 3) %>%
                mutate(
                    cutoff_y = cutoff,
                    model = "sharp",
                    covs = 0,
                    outcome = o,
                    men = ss_labs[i]
                )
        }) %>% reduce(rbind)
    }) %>% reduce(rbind)
}) %>%
    reduce(rbind) %>%
    mutate(
        p_order = paste0("Polynomial Order = ", p),
        method = "rbc",
        bw_select = "auto",
        east = 1
    )

## Naming

outcomevars <- c(
    "hh_income", "unemployed"
)

labels <- c(
    "Household income", "Unemployed"
)

## ##

res_main_m <- res_main_m %>%
    filter(outcome %in% c(outcomevars)) %>%
    left_join(data.frame(
        outcome = outcomevars,
        outcome_lab = labels,
        stringsAsFactors = F
    ))

## rename cohort cutoff

res_main_m <- res_main_m %>%
    mutate(cutoff_y = paste0("Cohort cutoff: ", cutoff_y))

## Plot this

pd <- position_dodge(0.5)

# Plot

vars_use <- c("Household income", "Unemployed")

p_df <- res_main_m %>%
    filter(str_detect(cutoff_y, "1971")) %>%
    filter(outcome_lab %in% vars_use) %>%
    mutate(men = ifelse(men == "Women", "East German\nwomen",
        "East German\nmen"
    ))

## Plot this

pd <- position_dodge(0.5)

p1 <- ggplot(
    p_df,
    aes(outcome_lab, estimate, covs)
) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "grey60") +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
        position = pd, width = 0
    ) +
    geom_point(
        shape = 21, fill = "white",
        position = pd, size = 2
    ) +
    theme_bw(base_size = 16) +
    coord_flip() +
    scale_color_manual(values = c("black", "grey40"), name = "") +
    scale_fill_manual(values = c("black", "grey40"), name = "") +
    scale_shape_manual(values = 21:22, name = "") +
    theme(legend.position = "bottom") +
    ylab("Effect of one additional year of\neducation under autocracy") +
    xlab("") +
    facet_wrap(~men)
p1
